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I. INTRODUCTION 

With a wide variety of practical applications and numerous opportunities for investigating 
basic questions in chemistry and physics, colloidal systems (suspensions of solid particles) 
and emulsions (dispersions of surfactant-coated liquid droplets) are of great interest to a 
variety of researchers. Recently, there has been increased attention given to emulsions 
involving water droplets dispersed in liquid crystals Surfactant-coated water droplets 
dispersed in a nematic host (so-called inverted emulsions, as opposed to direct emulsions 
where the nematic itself is separated into spherical droplets) can serve as generators of 
topological defects because of the anchoring of the liquid crystalline molecules to the droplet 
surfaces. Sufficiently strong homeotropic anchoring at the surface of the droplet induces a 
radial hedgehog defect centered on the droplet. To compensate for this induced defect, 
another defect must be created in the medium. 

There are two defect configurations which are theorized to arise in this case — the hyper- 
bolic dipole and the Saturn ring (see Fig. |l|). In the dipole case, the radial hedgehog 
defect is compensated by another point defect (a hyperbolic hedgehog); the two together 
form the dipole. The Saturn ring is a hyperbolic wedge disclination line which encircles the 
droplet equatorially. At large distances, the Saturn ring has the appearance of a hyperbolic 
hedgehog. In both cases then, the radial hedgehog with defect index +1 is compensated by 
a hyperbolic defect with index —1. Note that in the dipole case, the line joining the droplet 
and the hyperbolic monopole is parallel to the director fi, while in the Saturn ring case, 
the normal to the plane of the ring is parallel to n. Theoretical and numerical work based 
on elastic theory |0J|,|,0 suggests that the dipole configuration is stable for micron sized 
droplets and sufficiently strong surface anchoring, while a Saturn ring should appear if the 
droplet size is reduced. If the strength of the surface anchoring is reduced a third director 
configuration can arise, termed the surface-ring configuration [^-0. In this configuration 
the molecules at the equator are actually parallel to the droplet 's surface (homotropic) , and 
radial anchoring is limited to the poles of the droplet. Experiments to-date have always 



I INTRODUCTION 



3 



revealed dipolar configurations [|l|, consistent with the relatively large particle sizes and 
surface anchoring used. Changing the surfactant compound to produce homotropic (tan- 
gential) anchoring at the surface results in a pair of point defects at the poles of the droplet 
(boo j urns). 

The dynamics of the dispersed droplets and their accompanying director configurations 
are also of interest @,|||. Detailed knowledge of the fiow properties of liquid crystals is 
important for both fundamental and technological reasons. One way to measure viscosities 
in a fiuid is via the "falling ball experiment" . The drag force on a falling ball can be related 
to the viscosity by a Stokes-type formula: = —GtcRvi], where R is the radius of the 
ball, V is its terminal velocity and 1] is the fiuid viscosity. Because of the anisotropy of liquid 
crystalline molecules the drag force on a moving water droplet in an inverse nematic emulsion 
will depend on the relative orientation of the velocity and the director, i.e., a Stokes-like 
formula incorporates a nematic resistance tensor M 

F = Mxv + (M|, -Mx)(v -11)11, (1) 

where the subscripts _L and || refer respectively to directions perpendicular and parallel to the 
director. Eqn (|I]) implies that the drag force will not be antiparallel to the velocity, but will 
instead have a component, Fuft perpendicular to v. The drag force is particularly sensitive 
to the defect structure around the droplet because regions of high gradients in the director 
field (such as those around defects) result in higher resistance to fiow. Theoretical work on 
the drag force exerted on spheres as well as cylinders has been carried out by Ruhwandl and 
Terentjev in the low Ericksen number regime (small values of vRrj/K where v and R 
are the droplet velocity and radius respectively and r] and K are characteristic values of the 
nematic viscosity and elastic constant respectively) where they found a drag force anisotropy 
(the ratio Fj_/F||) in the range of 1.5-1.75. A drag force anisotropy not equal to unity is 
equivalent to a nonzero lift force. These authors also considered discotic molecules 0] and 
found a ratio less than unity which is reasonable given the oblate shape of the molecule. 
Here we present the results of molecular dynamics (MD) simulations carried out on a 
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system consisting of one spherical droplet immersed in a Gay-Berne nematic liquid crystal 
The Gay-Berne model is a phenomenological single site model of mesogenic molecules 
which captures the essential physical features of a wide variety of liquid crystalline behavior 
p!0| . Using a fluid model like the Gay-Berne system, as opposed to a lattice model of liquid 
crystals, we are able to simulate not only static behavior such as the director configuration 
surrounding the sphere (including the role of thermal fluctuations which are neglected in 
elastic theory) but dynamical behavior as well. We model the interaction of the Gay-Berne 
molecules with the droplet using a generalized repulsive Gay-Berne potential Jill account 



for excluded volume effects along with a phenomenological surface anchoring potential. Due 
to the large computational overhead of the Gay-Berne system our study is restricted to 2048 
Gay-Berne particles and a relatively small droplet, with a diameter equal to the long axis of 
the Gay-Berne molecule. However, using massively parallel computers it is possible to study 
much larger Gay-Berne systems consisting of at least 65,000 molecules |1l2| , p!3| and we hope 
to carry out a simulation of the emulsion system with a similar number of host molecules in 
the future. 

The outline of this paper is as follows. In the next section we present the computational 
details of our study. Section III contains the results for the director configuration in the 
absence of an external force on the droplet, as well as the results of "falling ball" type 
experiments. We offer some concluding remarks in Section IV. 

II. COMPUTATIONAL DETAILS 



We model the nematic host using the original Gay-Berne potential given by ||14|| : 

U (ui,Uj,f) = 4e(ui,Uj,f) 



>, 12 f 

O"0 O-Q 



(2) 



r - a (ui, Uj, f) + (To J [r - a (ui, uj, f) + ao J 
where Ui, uj give the orientations of the long axes of molecules i and j, respectively, and r 
is the intermolecular vector (r = ri — rj). The parameter a (ui, uj, f) is the intermolecular 
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separation at which the potential vanishes, and thus represents the shape of the molecules. 
Its explicit form is 

\2 



cr (Ui,Uj,f) = (To 



r ■ Ui + r ■ u; 



2-" [ 1 + X(ui-Uj) 

2 ^ 1 -1/2 



+ 



where ctq = (defined below) and x is 



r • U; — r ■ U; 



1 - X (ui ■ Uj 



(3) 



(4) 



Here is the separation between two molecules when they are oriented end-to-end, and ag 
the separation when side- by-side. The well depth e (ui, uj, f), representing the anisotropy of 
the attractive interactions, is written as 



s (ui, Uj, f) = eoe'' (u;, uj) e'^" (ui, uj, f) 



(5) 



where 



£(Ui,Uj) = {l-x^ (Ui -Uj)^} 



'1/2 



(6) 



and 



fr- - -\ 1 1 / [ (f ■ Ui + r ■ Uj) 



+ 



r ■ Ui - r ■ u. 



(7) 



l-X'(ui-Uj) J 

with x' defined in terms of Se and Eg, the end-to-end and side-by-side well depths, respec 
tively, as 



X' = {l-{ee/esf'}/{l + {ee/esf^} 



(8) 



The overall energy scale is set by the value of Eq- 

For the adjustable parameters, we primarily used the values originally proposed by Gay 
and Berne @: n = 2,1/ = X^Of^jos = 3, and Es/e^ = 5. We also ran simulations with the 
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parameterization suggested by Berardi et al. (where fi = 1 and z/ = 3) and saw no 
noticeable differences in our results. 

The interaction between the Gay-Berne molecules and a droplet of radius R consists 
of two parts — a soft, repulsive term Ur and an anchoring term Ua- The repulsive term 
utilizes the parameters of a Gay-Berne potential generalized to mimic the interaction between 
nonequivalent particles Jill. The range parameter ar (u, f) and the energy parameter Er (u) 



used in this potential are generalizations of the corresponding parameters in Eqns. (^) and 
(I), namely 



a,(u,f) = <(l-X.(f -u))-'/' (9) 



and 



£,(u,f) = l, (10) 



with Xr defined as 



and 



^''J+^R^' ^^^^ 



^ (af + 4i?2)i/2. (12) 



Thus the repulsive part of the potential between a molecule with orientation u and the 
droplet with center of mass separation r is given by 

[/,(u,f)=4eoX . (13) 

[r-ar (u, r) + (To J 

We use an exponent of 18 instead of 12 in this repulsive term to make the interaction 
"harder" . 

The short-ranged surface anchoring term is given by: 

Ua = -W^-^, (14) 



Ill RESULTS 



7 



where is a phenomenological anchoring coefficient. Only those molecules with r < R + cTg 
are considered to interact with the droplet via this surface anchoring. A positive value of 
W yields homeotropic anchoring, while a negative value of W gives tangential anchoring. 

We performed MD simulations at constant temperature, T* = ksT/eo = 0.95, and pres- 
sure, P* = Pa^/So = 4.0, corresponding to a nematic phase. We controlled temperature 
and pressure using the equations of motion of refs. [l^ and Jl^, which are modifications of 



the original Nose-Hoover approach [|I^ . The orientational degrees of freedom were parame- 
terized using quaternions. We used a variation [T^ of the leapfrog algorithm to solve these 
equations for a system consisting of 2048 Gay-Berne molecules and one spherical droplet 
with a timestep At* = 0.001, At* = (mal/eo)~^^'^, where m is the mass of the Gay-Berne 
molecule. The droplet diameter was chosen to be Sctq which is the length of the long axis of 
the Gay-Berne molecules in the host, and the mass of the droplet was chosen to be 100m, 
so that in the absence of an external force the sphere moves very slowly compared to the 
molecules. The surface anchoring parameter W was chosen to be 4000 (in dimensionless 
units) for weak anchoring and 35000 for strong anchoring. In terms of the dimensionless 
quantity C, = K/WR which measures the elastic deformation energy relative to the surface 
anchoring energy the former value of W corresponds to ^ ~ 0.4, while the latter cor- 
responds to a value of ^ ten times smaller (these values assume that K is of order unity). 
We prepared our system by placing the sphere at the center of an empty MD cell. We 
then surrounded the sphere with eight replicas of a Gay-Berne nematic phase of 256 par- 
ticles, equilibrated to the temperature and pressure specified above. The system was then 
equilibrated for 100000 timesteps. 



III. RESULTS 
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A. Static Director Configuration 



In the absence of an external force on the droplet we studied the director configuration 
using a variety of techniques and considered both strong and weak surface anchoring. Af- 
ter the equilibration process described in the previous section we carried out production 
runs of 60000 steps examining the director configuration every 6000 steps. The general fea- 
tures of the configurations required to assess the overall director structure were essentially 
time-independent. We studied these configurations using the defect-finding algorithm for 
disclinations described in [|2T|, and contour plots of the angle of molecular orientation with 
respect to the director. To apply the defect-finding algorithm we divided the system into a 
6x6x6 lattice of cubic bins (the size is chosen so that the defects ultimately found have a 
core size of one bin). Note that the creation of a lattice is strictly for convenience in defect 
finding; the time evolution of the system allows for complete translational freedom. Within 
each bin, the order parameter tensor 



was calculated, its largest eigenvalue taken as the local order parameter S and the corre- 
sponding eigenvalue as the local director n. Topological defects can be located by considering 
the directors at the corners of a square (one of the faces of a cube) in the 3D lattice. The 
idea is to track the course of the intersections of these directors with the order parameter 
sphere (actually the projective plane RP2) as one moves around the corners of the real-space 
square (Fig. |^). Starting with the intersection of with the sphere, one then takes as the 
next point either the intersection of or — n^, whichever is closest to fill's intersection. 
Once this point is determined, either nc or —fie is used, depending on the proximity to the 
previously defined point, and so on. Once the last point (from corner D) is determined, one 
looks at whether its intersection is in the same hemisphere as the starting point's. If so, no 
defect is present — the path in order parameter space is deformable to a single point, i.e., 
a uniform configuration. If the first and last points are in different hemispheres, however, 




(15) 
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then a disclination line is taken to cut through the center of the square and is oriented 
perpendicular to the plane of the square. 

For strong anchoring [W = 35000,^ = 0.04) we detected a Saturn ring configuration 
using this algorithm (Fig. |^). Given the small size of our droplet this result is consistent 
with the results of elastic theory The molecular configuration in the portion of 

the MD cell surrounding the droplet is shown in Fig. § providing a clear view of the strong 
radial anchoring of the molecules immediately surrounding the sphere. It is clear that the 
radial anchoring takes place over one molecular layer. Thus, it is not possible with our 
system size to provide more than an estimate of the radius of the Saturn ring, namely, it 
is slightly greater than the radius of the droplet. However, our result is consistent with the 
results of elastic theory analyses of references ||^ and [^, who find radii of order R. With 
a larger system size and thus a larger droplet, it should be possible to determine the radius 
of the ring. 

The director configurations can also be visualized using contour plots of the angle of 
orientation with respect to the director. As shown in Fig. ^ (a slice though the center of the 
droplet and perpendicular to n), the molecular orientation is along the director everywhere 
except in a ring about the droplet. Fig. ^ shows a slice parallel to fi. Both figures are 
consistent with the presence of a Saturn ring in the equatorial plane perpendicular to n. A 
similar set of figures for tangential anchoring (Figs. ^ and show a pair of boojums located 
near the poles of the droplet. 

For the lower value of the anchoring parameter that we studied, W = 4000, we found 
that the Saturn ring is replaced by the surface-ring configuration as shown in Fig. ^. This 
value of W corresponds to a surface extrapolation length ^ = 0.4, and our result is consistent 
with the predictions of elastic theory [p|,|22|. 
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B. Dynamics 

To study the dynamics of the immersed droplet we began with an equihbrated config- 
uration prepared as described at the end of section II. We applied a force to the droplet 
and measured its position as a function of time, determining the terminal velocity from the 
asymptotic slope of the droplet's position as a function of time (see Fig. 0). For a force 
applied at an angle of 45° with the director the lift force ratio (the ratio of the components 
of drag force perpendicular and parallel to the velocity) is given by, 

= tan(45°-^), (16) 

where 6 is the angle between the terminal velocity vector and the director. The anisotropy 
of the resistance tensor M is given by, 

M^/Mii=vii/v^ = cote. (17) 

Table | shows the magnitude of the terminal velocity, the value of 6, the lift force ratio and the 
resistance tensor anisotropy for driving forces of different magnitudes for both W = 35000 
(strong anchoring) and for W = 0. The magnitudes of M± and My in the former case are 
approximately five times greater than in the latter, as one might expect intuitively. For both 
cases, the values of 6 are generally slightly less than 45°. Thus, the lift force ratio is nearly 
zero, and the resistance anisotropy nearly unity. These results are very different from what 
is found experimentally and predicted theoretically for droplets which are many times 
larger than the nematic host molecules (the theoretical work uses the Frank elastic energy, 
and thus assumes large droplets). However, given our small droplet size (diameter equal to 
the length of a nematic molecule), it is not surprising that the lift force Fuji is very small, 
and the resistance tensor nearly isotropic. 

In the anchoring case, there is very little distortion of the Saturn rings for small driving 
forces. With larger forces, however, the Saturn ring either folds back on itself behind the 



droplet for parallel driving forces (Figs. ITT] and O) or becomes stretched out behind the 
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droplet while retaining its equatorial position for perpendicular driving forces (Figs. |T3| and 
[M]). The distortions become so large that, as one can clearly see in Fig. |12|, for example, 
when the droplet exits the left side of the simulation box and re-enters on the right, it 
will begin to interact with its own wake. For the no-anchoring case, there is no Saturn 
ring to distort, but the droplet still causes localized disturbances which eventually result in 
self-interactions for large driving forces. 

IV. CONCLUSIONS 

We have shown in an off-lattice simulation that, consistent with theoretical predictions, 
the Saturn ring is the preferred defect configuration around small droplets with radial an- 
choring embedded in a nematic. We have also shown that lowering the anchoring strength 
can produce a surface ring instead of a Saturn ring. The Saturn ring increases (compared 
to the no-anchoring case) resistance to the droplet's motion when the droplet is subjected 
to a driving force. We found a very small lift force (and correspondingly little anisotropy in 
the resistance tensor), consistent with the small droplet size studied. Clearly, further work 
with larger system sizes will allow numerous interesting questions to be studied fruitfully. 
The size limitation on the droplet comes, of course, from periodic boundary conditions ; if 
the droplet is too large, the anchoring layers on opposite sides of the droplet will begin to 
interact with each other. Larger system sizes will also allow exploration of the variation in 
Saturn ring radius with particle size and, very importantly, study of the transition from Sat- 
urn ring to dipole. With larger systems, one can also reasonably begin to study multiparticle 
dispersions. 
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FIGURES 

FIG. 1. Two possible defect configurations surrounding a water droplet placed in a nematic 
host, with strong radial surface anchoring. The droplet is indicated by the large sphere, (a) Dipole 
configuration showing the induced hyperbolic monopole, indicated by the small sphere, (b) Saturn 
ring consisting of an induced wedge disclination ring about the droplet. The Saturn ring extends 
into and out of the plane of the page. 

FIG. 2. The disclination- finding algorithm. The directors at the corners of the lattice cube 
face shown on the left are tracked on the order parameter space sphere shown on the right. The 
diameters AA' , BB' ,CC' , DD' correspond to the axes of the headless director at the real-space 
lattice sites A, B, C, D respectively. 

FIG. 3. A Saturn ring surrounding the droplet (of diameter Sctq) at the center of the MD cell, 
with surface anchoring parameter W = 35000 (strong anchoring). The ring was located using 
the disclination-finding algorithm described in The MD cell was divided into a 6 x 6 x 6 

lattice leading to the jagged appearance of the ring. Most of the ring segments are located in the 
equatorial plane perpendicular to the director. 

FIG. 4. Molecular configuration for strong anchoring. A subset of the total simulation box 
centered on the droplet is shown and those molecules closest to the droplet are highlighted to 
assist in viewing the defect configuration. 

FIG. 5. Contour plot of the director configuration in a slice centered on a droplet with strong 
radial anchoring. The lighter regions represent molecules aligned with the director, while darker 
regions indicate molecules perpendicular to the director. Here the director points up out of the 
page, so the figure shows a Saturn ring in the plane of the figure. 

FIG. 6. Contour plot of the director configuration for strong radial anchoring in a slice with 
the director in the plane of the page as indicated. The lighter regions represent molecules aligned 
with the director, while darker regions indicate molecules perpendicular to the director. 
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FIG. 7. Contour plot of the director configuration with parallel surface anchoring in a slice 
perpendicular to the director. Compare with the case of radial anchoring shown in Fig. |5|. 

FIG. 8. Contour plot of the director configuration with parallel surface anchoring in a slice in 
the plane of the director. Compare with the case of radial anchoring shown in Fig. ^. 

FIG. 9. The same as Fig. ^ but with weaker radial anchoring {W = 4000). The Saturn ring 
has contracted to a surface ring. 

FIG. 10. Illustration of asymptotic velocity measurement. Sphere positions x,y,z (top to 
bottom curves respectively) are shown for a driving force F = 85 applied at 45° to the director. 
The slopes of the straight lines yield the corresponding asymptotic velocity components Vx,Vy,Vz- 

FIG. 11. Slice perpendicular to the director for driving force F = 250, with the force applied 
parallel to the director. 

FIG. 12. Slice in the plane of the director for driving force F = 250, with the force applied 
parallel to the director. The Saturn ring is clearly distorted. 

FIG. 13. Slice perpendicular to the director for driving force F = 250, with the force applied 
perpendicular to the director. Again, the Saturn ring is clearly distorted. 

FIG. 14. Slice in the plane of the director for driving force F = 250, with the force applied 
perpendicular to the director. 
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TABLES 

TABLE I. Lift force ratio and resistance tensor anisotropy. A force of magnitude F is applied 

at an angle of 45° to the director and the magnitude v and angle 6 (with respect to the director) 
of the terminal velocity are measured. The lift force ratio Fnft/Fy is given by tan(45° — 6) and the 
resistance tensor anisotropy M±/M\\ by cot^. . 
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